function [ mu, v ] = biv_logn( mean_mu, var_mu, mean_v, var_v, rho, K)


% produces sets of consumers for arbitrary market
% the expressions for the conditional distribution of mu can be found,
% https://en.wikipedia.org/wiki/Multivariate_normal_distribution
% the inputs are the parameters of the corresponding GAUSSIAN
% types are computed by quadrature using the inverse CDF of the lognormal - ie, types are not random
% the input is K, but K^2 individuals are simulated
    
%% safeguard

assert(-1<=rho && rho<=1, 'must have -1<=rho<=1')

% if rho>1
%     error(['rho=' ns(rho) '>1']); 
% end

%% market statistics

sd_v = sqrt(var_v);

%% generating consumers

mu  = NaN(K^2,1);                                   % pre-allocating memory for K^2 consumers with 2D types
v   = NaN(K^2,1);
uniform_points = linspace( 1/K, 1-(1/K), K)';                   % vector of K evenly spaced points within (0,1)
v_points = logninv(uniform_points, mean_v, sd_v);                 % baseline vector of v's, will be repeated K times
for k=1:K       % for each value of v, generate corresponding values of mu conditional on v 
    v_k = v_points(k);
    mean_mu_given_v      = mean_mu + rho * sqrt(var_mu/var_v) * ( log( v_k ) - mean_v ) ;     % E[mu| v]
    var_mu_given_v       = (1 - rho^2) * var_mu;                                                 % V[mu|v]
    sd_mu_given_v        = sqrt(var_mu_given_v);
    mu_k = logninv( uniform_points, mean_mu_given_v, sd_mu_given_v );      % vector of mu, lognormal with mean E[mu|v] and corresponding standard deviation
    v(  K*(k-1)+1 : K*k, 1) = v_k*ones(K,1) ;
    mu( K*(k-1)+1 : K*k, 1) = mu_k ;                   % stacking up the results to obtain K^2 consumers                                                  % for each value of mu_k, K values of v are generated                                         
end


end

